The objectives of the exercises of this assignment are:
1) Download and organise the data and model and plot staircase responses based on fits of logistic functions
2) Fit multilevel models for response times
3) Fit multilevel models for count data
REMEMBER: In your report, make sure to include code that can reproduce the answers requested in the exercises below (MAKE A KNITTED VERSION)
REMEMBER: This assignment will be part of your final portfolio
Go to https://osf.io/ecxsj/files/ and download the files associated with Experiment 2 (there should be 29).
The data is associated with Experiment 2 of the article at the following DOI https://doi.org/10.1016/j.concog.2019.03.007
datadir<-"/Users/Laura/Desktop/Methods 3/github_methods_3_class/week_03/experiment_2"
laurafiles<-list.files(datadir,pattern='.csv',full.names=TRUE)
df <- vroom(laurafiles)
## Rows: 18131 Columns: 17
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): trial.type, task, target.type, obj.resp, subject
## dbl (12): pas, trial, jitter.x, jitter.y, odd.digit, target.contrast, target...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
obj.resp [i,16], o, e target_type [i,11], odd, even correct [i,19]
df$correct <- ifelse((df$obj.resp == "o" & df$target.type=="odd") | (df$obj.resp == "e" & df$target.type=="even"), 1,0)
class they should be classified into, e.g. factor, numeric etc.trial.type: a character variable encoded as either “staircase” (i.e. data was collected as a part of the adaptive staircase procedure which was performed at the beginning of the study) or “experiment” (i.e. data was collected as a part of the actual experiment trials). It should be re-encoded a factor in order for trial.type to be modelled later in the analysis, if relevant. However, if we were only interested in the experimental trials, we might consider filtering out all staircase trials*.
df$trial.type <- as.factor(df$trial.type)
pas: indicates one of the 4 categorically different ratings of the Perceptual Awareness Scale (PAS: Ramsøy & Overgaard, 2004), rated by the participant. It indicates the participant’s perceptual awareness of the stimulus where the numbers represent the following:*
(1) No experience (i.e. "no impression of the simulus. All answes are seen as mere guesses)*
(2) Weak Glimpse (i.e. "A feeling that something has been shown. Not characterized by any content, and this cannot be specified any further)*
(3) Almost Clear Experience (i.e. “Ambiguous experience of the stimulus. Some stimulus aspects are experienced more vividly than others. A feeling of almost being certain about one’s answer”)*
4) Clear experience (i.e. “Non-ambiguous experience of the stimulus. No doubt in one’s answer”)*
It should be encoded as a factor as it is a case of ordinal data.
df$pas <- as.factor(df$pas)
trial: refers to the trial number. It is numeric but should be factor
df$trial <- as.factor(df$trial)
target.contrast: indicates the contrast of the target stimulus, relative to the background, adjusted to match the threshold of each individual by using the QUEST-algorithm. It is numeric and should stay numeric as we are dealing with a continuous variable with values falling within the range of 0-1*
cue: cue shown to participants to indicate the set of possible digits that might appear on each trail. There were 36 different combination of cues (0-35). The task-column specifies the number of potential targets cued. As the cue variable contains nominal data, it should be reencoded into a factor.
df$cue <- as.factor(df$cue)
task: indicates the number of potential targets that were cued. Consisted of either 2,4, or 8 digits, indicated by the character variable task taking one of the 3 values:
singles (2 possible targets) e.g. 2:9 pairs (4 possible targets) e.g. 24:57 quadruplet (8 possible targets) e.g. 2468:3579
This variable should also be re-encoded as a factor as we are, again, dealing with nominal data.
df$task <- as.factor(df$task)
target_type: indicates the parity of the target stimulus (i.e. whether the digit was even or odd). It is encoded as a character but should be re-encoded as a factor.
df$target.type <- as.factor(df$target.type)
rt.subj: response time of answering Perceptual Awareness Scale Encoded as numeric as it is continuous data.
rt.obj: response time (time taken by participant to indicate the parity of the target stimulus). Encoded as numeric as it is continuous data.
obj.resp: the participant’s answer when asked to indicate the parity of the target stimulus (e = even, o = odd). Encoded as numeric but should be re-encoded as a factor.
df$obj.resp <- as.factor(df$obj.resp)
subject: participant ID, enabling us to distinguish between data collected from different participants. Is a character variable but should be re-encoded as a factor as we are dealing with a repeated measures design and therefore need to be able to model subject as a random effect in the analysis phase.
df$subject <- as.factor(df$subject)
correct: indicates accuracy of the respondent’s answer of the target stimulus’ parity. It is a binary variable where 0 indicates an incorrect response and 1 indicates a correct one, and should therefore be re-encoded as a factor.
df$correct <- as.factor(df$correct)
glm) that models correct as dependent on target.contrast. These plots will be our no-pooling model. Comment on the fits - do we have enough data to plot the logistic functions?library(tidyverse)
# make subset of df with only trial.type == staircase
df_staircase <- subset(df, df$trial.type == "staircase")
df_staircase$subject = gsub("(?<![0-9])0+", "", df_staircase$subject, perl = TRUE)
df_staircase$subject = as.integer(df_staircase$subject)
nopoolfun <- function(i){
dat <- df_staircase[which(df_staircase$subject == i),]
model <- glm(correct~target.contrast, family = 'binomial', data=dat)
fitted <- model$fitted.values
plot_dat <- data.frame(cbind(fitted,'target.contrast'=dat$target.contrast))
plot <- ggplot(plot_dat, aes(x = target.contrast, y = fitted))+
geom_point()+
geom_line(aes(x = target.contrast, y = fitted))+
xlab('Target Contrast')+
ylim(c(0,1))+
theme_minimal()
print(plot)
}
# Running the function for each participant
for (i in 1:29){
nopoolfun(i)
}
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
subjects <- c(1:16)
plots <- lapply(subjects, FUN=nopoolfun)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
do.call(grid.arrange, plots)
subjects <- c(17:29)
plots <- lapply(subjects, FUN=nopoolfun)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
do.call(grid.arrange, plots)
# --------------- No pooling model ---------------
m_nopool <- glm(correct~target.contrast*subject, family=binomial(link=logit), data=df_staircase)
fitted_nopool <- fitted(m_nopool)
df_staircase$fitted_nopool <- fitted_nopool
# one plot
ggplot(data = df_staircase)+
geom_point(aes(x = target.contrast, y = fitted(m_nopool), color = subject))+
geom_line(aes(target.contrast, fitted(m_nopool))) +
labs(title = "No pooling plot")+
facet_wrap(.~subject)
From visual inspection of the no pooling plot, it becomes clear that the fits are very bad. A subset data frame containing only the staircasing data doesn’t contain enough data to plot the logistic functions (the fitted values do not follow a sigmoid function).
glmer from the package lme4) where unique intercepts and slopes for target.contrast are modelled for each subject# Telling the model to expect differing baseline-levels of "correct" (the intercept, represented by 1) as well as differing responses to the main factor in question, which is "target.constrast" in this case.
m_partial <- glmer(correct~target.contrast + (1+target.contrast|subject), data = df_staircase, family = "binomial")
fitted_partial <- fitted(m_partial)
df_staircase$fitted_partial <- fitted_partial
## Plotting on top
ggplot(data = df_staircase)+
geom_point(aes(x = target.contrast, y = fitted(m_nopool), color = "no pooling"))+
geom_point(aes(x = target.contrast, y = fitted(m_partial), color= "partial pooling"))+
geom_line(aes(target.contrast, fitted(m_partial), color = "partial pooling")) +
geom_line(aes(target.contrast, fitted(m_nopool), color = "no pooling")) +
labs(title = "No pooling and partial pooling model plot")+
facet_wrap(.~subject)
Partial pooling allows for a better fit for each subject as it takes into account the expected individual baseline differences depending on the stimulus contrast level as well as allows for different slopes for each subject. In other words, it allows the model to reflect that individuals would perform differently in different experimental settings. Partial pooling provides us a compromise between complete and no-pooling models, as it allows us to model both an average (general tendency of the whole dataset) and each level of the categorical predictor, subject. In partial pooling (a mixed effects model), we model both an average and each level, i.e. combine information from the population (fixed) effects of the complete pooling model and the subject-specific (random) effects of the no-pooling one. (((accounts for nested structure of our data, it assumes hierarchy in our data structure)))
Now we only look at the experiment trials (trial.type)
df_experiment <- df %>% filter(trial.type=="experiment")
df_experiment
## # A tibble: 12,528 x 18
## trial.type pas trial jitter.x jitter.y odd.digit target.contrast
## <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 experiment 1 0 -0.322 0.395 5 0.0392
## 2 experiment 1 1 0.499 0.495 3 0.0392
## 3 experiment 1 2 0.461 -0.198 5 0.0392
## 4 experiment 1 3 0.200 -0.0263 3 0.0392
## 5 experiment 1 4 -0.347 -0.0448 5 0.0392
## 6 experiment 1 5 0.0944 -0.106 5 0.0392
## 7 experiment 1 6 0.291 -0.347 5 0.0392
## 8 experiment 1 7 0.00541 -0.168 3 0.0392
## 9 experiment 2 8 -0.393 -0.357 5 0.0392
## 10 experiment 1 9 -0.110 0.185 3 0.0392
## # … with 12,518 more rows, and 11 more variables: target.frames <dbl>,
## # cue <fct>, task <fct>, target.type <fct>, rt.subj <dbl>, rt.obj <dbl>,
## # even.digit <dbl>, seed <dbl>, obj.resp <fct>, subject <fct>, correct <fct>
df_experiment$subject = gsub("(?<![0-9])0+", "", df_experiment$subject, perl = TRUE)
df_experiment$subject = as.integer(df_experiment$subject)
# pick 4 random from data frame
library(dplyr)
set.seed(1)
sample <- sample_n(df_experiment,4) # note that you could get the same number twice
sample # 03, 019, 012, 025
## # A tibble: 4 x 18
## trial.type pas trial jitter.x jitter.y odd.digit target.contrast
## <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 experiment 1 152 0.362 0.479 7 0.0628
## 2 experiment 1 227 0.450 0.188 5 0.0357
## 3 experiment 1 22 0.161 -0.269 9 0.106
## 4 experiment 4 0 0.399 0.0117 5 0.0309
## # … with 11 more variables: target.frames <dbl>, cue <fct>, task <fct>,
## # target.type <fct>, rt.subj <dbl>, rt.obj <dbl>, even.digit <dbl>,
## # seed <dbl>, obj.resp <fct>, subject <int>, correct <fct>
# make df subsets for each subject
sub3 <- df_experiment %>% filter(subject==3)
sub19 <- df_experiment %>% filter(subject==19)
sub12 <- df_experiment %>% filter(subject==12)
sub25 <- df_experiment %>% filter(subject==25)
# make modes (lm cause continuous variables)
exp_model <- lm(rt.obj~1, data = df_experiment) # overall
# FModels
m3 <- lm(rt.obj~1, data = sub3)
m19 <- lm(rt.obj~1, data = sub19)
m12 <- lm(rt.obj~1, data = sub12)
m25 <- lm(rt.obj~1, data = sub25)
# QQ-plot
qq3 <- ggplot(sub3, aes(sample = residuals(m3))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subject 003: Q-Q Plot for residuals") +
theme_bw()
qq19 <- ggplot(sub19, aes(sample = residuals(m19))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subject 019: Q-Q Plot for residuals") +
theme_bw()
qq12 <- ggplot(sub12, aes(sample = residuals(m12))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subject 012: Q-Q Plot for residuals") +
theme_bw()
qq25 <- ggplot(sub25, aes(sample = residuals(m25))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subject 025: Q-Q Plot for residuals") +
theme_bw()
ggpubr::ggarrange(qq3, qq19, qq12, qq25)
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
i. comment on these
All 4 qq-plots indicates major deviances from normality, i.e. right skewness. We could try to improve the QQ-plots via log-transformation of the response time data.
m3log <- lm(log(rt.obj)~1, data = sub3)
qq3log <- ggplot(sub3, aes(sample = residuals(m3log))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subj 3: QQ-plot residuals (log)") +
theme_bw()
m19log <- lm(log(rt.obj)~1, data = sub19)
qq19log <- ggplot(sub19, aes(sample = residuals(m19log))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subj 19: QQ-plot residuals") +
theme_bw()
m12log <- lm(log(rt.obj)~1, data = sub12)
qq12log <- ggplot(sub12, aes(sample = residuals(m12log))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subj 12:QQ-plot residuals") +
theme_bw()
m25log <- lm(log(rt.obj)~1, data = sub25)
qq25log <- ggplot(sub25, aes(sample = residuals(m25log))) + stat_qq() +
stat_qq_line(colour = "red") +
labs(x = "Theoretical quantiles", y = "Sample quantiles") + ggtitle("Subj 25: QQ-plot residuals") +
theme_bw()
ggpubr::ggarrange(qq3log, qq19log, qq12log, qq25log)
Doing a log-transformation improves the QQ-plots slightly, altough the data still deviates from normality.
REML=FALSE in your lmer-specification)partial1 <- lmerTest::lmer(log(rt.obj)~task + (1|subject) + (1|trial), data = df_experiment, REML=FALSE)
partial2 <- lmerTest::lmer(log(rt.obj)~task + (1|trial), data = df_experiment, REML=FALSE)
partial3 <- lmerTest::lmer(log(rt.obj)~task + (1|subject), data = df_experiment, REML=FALSE)
partial4 <- lmerTest::lmer(log(rt.obj) ~ task + (1 + task | subject) + (1 | trial), data = df_experiment, REML = F)
partial5 <- lmerTest::lmer(log(rt.obj) ~ task + (1 + task | subject), data = df_experiment, REML=FALSE)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -7.6e+02
partial_model <- c("m1", "m2", "m3", "m4", "m5")
AIC <- c(AIC(partial1), AIC(partial2), AIC(partial3), AIC(partial4), AIC(partial5))
sigma <- c(sigma(partial1),sigma(partial2),sigma(partial3), sigma(partial4), sigma(partial5))
# calculating pseudo R-squared for mixed effects model (R2m = marginal R2, provides the variance explained only by fixed effects. R2c = conditional R2, provides the variance explained by the entire model, i.e., both fixed effects and random effects)
pacman::p_load(MuMIn)
r.squaredGLMM(partial1) # R2c = .2629
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
## R2m R2c
## [1,] 0.006697664 0.2628567
r.squaredGLMM(partial2) # R2c = .0360
## R2m R2c
## [1,] 0.006663833 0.03597408
r.squaredGLMM(partial3) # R2c = .2253
## R2m R2c
## [1,] 0.006538142 0.2253338
r.squaredGLMM(partial4) # R2c = .2750
## R2m R2c
## [1,] 0.006742469 0.27504
r.squaredGLMM(partial5) # R2c = .2371
## R2m R2c
## [1,] 0.006538065 0.2370859
pseudo_r2 <- c(0.2629, 0.0360, 0.2253, 0.2750, 0.2371)
as_tibble(cbind(partial_model, sigma, AIC, pseudo_r2))
## # A tibble: 5 x 4
## partial_model sigma AIC pseudo_r2
## <chr> <chr> <chr> <chr>
## 1 m1 0.767353686788801 29459.9382110624 0.2629
## 2 m2 0.877443929282292 32560.1503149819 0.036
## 3 m3 0.786510836280231 29685.3126188651 0.2253
## 4 m4 0.761030254763385 29327.6815811647 0.275
## 5 m5 0.780526731035549 29560.3359954543 0.2371
m4 has the lowest sigma as well as AIC value. This is the model that predicts objective response times dependent on task (fixed effects), modelling a random intercepts for subjects and trial as well as random slopes for task. Random intercepts were modeled for subject as one wouold expect individual baseline performance differences among different participants. Trial as also been modelled as random intercepts as you would expect different trials to vary in difficulty (and participants might modify their performance as they get acquainted with the experiment).
summary(partial4)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: log(rt.obj) ~ task + (1 + task | subject) + (1 | trial)
## Data: df_experiment
##
## AIC BIC logLik deviance df.resid
## 29327.7 29409.5 -14652.8 29305.7 12517
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.5363 -0.4958 -0.0279 0.5119 8.0799
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## trial (Intercept) 0.029952 0.17307
## subject (Intercept) 0.149162 0.38621
## taskquadruplet 0.004001 0.06325 0.35
## tasksingles 0.032270 0.17964 0.38 -0.58
## Residual 0.579167 0.76103
## Number of obs: 12528, groups: trial, 432; subject, 29
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.26407 0.07317 29.75753 -3.609 0.00111 **
## taskquadruplet -0.07225 0.02054 28.60783 -3.518 0.00147 **
## tasksingles -0.17868 0.03739 29.12028 -4.778 4.65e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tskqdr
## taskqudrplt 0.099
## tasksingles 0.279 -0.114
The estimates for quadruplet and singles are negative, which signifies that these tasks have a shorter reaction time compared to the pair task. Relative to the quadruplet the singles task is generally faster.
The coefficients “taskquadruplet” and “tasksingles” are both negative and accompanied by small significant p-values, indicating that as the task changes from pairs to taskquadruplet and Tasksingles, respectively, the response time significantly decreases. The effect is, however, strongest from changing from pairs to singles (the respodents generally complete the singles tasks faster, relative to the quadruplet tasks).
# fails to converge: pastask <- lmer(log(rt.obj)~pas*task + (1 + task | subject) + (1 | trial), data = df_experiment, REML=FALSE)
pastask <- lmer(log(rt.obj)~pas*task + (1 | subject) + (1 | trial), data = df_experiment, REML=FALSE)
pastask1 <- lmer(log(rt.obj) ~ task*pas + (1|subject), data = df_experiment, REML = FALSE)
pastask2 <- lmer(log(rt.obj) ~ task*pas + (1|trial) + (1|subject), data = df_experiment, REML = FALSE)
pastask3 <- lmer(log(rt.obj) ~ task*pas + (1|trial) + (1|subject) + (1|target.contrast), data = df_experiment, REML = FALSE)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
pastask4 <- lmer(log(rt.obj) ~ task*pas + (1|trial) + (1|subject) + (1|cue), data = df_experiment, REML = FALSE)
pastask5 <- lmer(log(rt.obj) ~ task*pas + (1|trial) + (1|subject) + (1|cue) + (1|target.contrast), data = df_experiment, REML = F)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
The first 3 models (up to 3 group intercept) can be modelled without running into convergence issues if I add trial, subject and cue. However, if I try to add target.contrast as my third group intercept, I run into convergence issues. However, if adding cue as my third group intercept and target.contrast at my fourth, I can add 3 types (but run into issues at the fourth one).
print(VarCorr(<your.model>), comp='Variance') to inspect the variance vector - explain why the fit is singular (Hint: read the first paragraph under details in the help for isSingular)print(VarCorr(pastask5), comp='Variance')
## Groups Name Variance
## trial (Intercept) 0.0252761
## cue (Intercept) 0.0025684
## target.contrast (Intercept) 0.0644419
## subject (Intercept) 0.0930045
## Residual 0.5592559
help(isSingular)
We can see that the variances of all the random intercepts in the model pastask5 are all very low (especially our random-effect variance estimate for “cue” is nearly zero). Variances being estimated as 0 results in singular fit. It indicates that the model is overfitting, so the structure of our random effects is too complex for our data to support it.
data.count. count should indicate the number of times they categorized their experience as pas 1-4 for each task. I.e. the data frame would have for subject 1: for task:singles, pas1 was used # times, pas2 was used # times, pas3 was used # times and pas4 was used # times. You would then do the same for task:pairs and task:quadrupletpacman::p_load(dplyr)
data_count <- df %>%
group_by(subject, task, pas) %>%
dplyr::summarize("count" = n())
## `summarise()` has grouped output by 'subject', 'task'. You can override using the `.groups` argument.
data_count
## # A tibble: 340 x 4
## # Groups: subject, task [87]
## subject task pas count
## <fct> <fct> <fct> <int>
## 1 001 pairs 1 109
## 2 001 pairs 2 45
## 3 001 pairs 3 4
## 4 001 pairs 4 12
## 5 001 quadruplet 1 141
## 6 001 quadruplet 2 23
## 7 001 quadruplet 3 3
## 8 001 quadruplet 4 1
## 9 001 singles 1 43
## 10 001 singles 2 87
## # … with 330 more rows
# the cause of this error is R’s confusion which summarize function (dplyr vs. plyr) it should use.
countmodel1 <- glmer(count~pas*task + (1+pas|subject), data = data_count, family = poisson, control = glmerControl(optimizer="bobyqa"))
summary(model1) i. which family should be used?
Poisson Regression should be used as it is best for modeling a dependent variable that consists of “count data” given one or more independent variable. Poisson distributed data is integer-valued, discrete, not continuous, and is limited to non-negative values. It also assumes that the errors follow a Poisson distribution rather than a normal distribution, and models the natural log of the response variable, ln(Y), as a linear function of the coefficients (rather than modeling Y as a linear function of the regression coefficients).
pas is encoded as a factor, meaning that it computes the analysis for each level seperately . Thus, we’ll get a slope for each separate level of pas, where each slope is relative to the reference level.
dfoptim package is needed. In glmer, you can add the following for the control argument: glmerControl(optimizer="bobyqa") (if you are interested, also have a look at the function allFit)The control argument “glmerControl(optimizer=”bobyqa“)” was added to avoid convergence error
# Model with only the main effects of pas and task
countmodel2 <- glmer(count ~ pas + task + (1 + pas|subject), data = data_count, family = "poisson", control = glmerControl(optimizer="bobyqa"))
rownames <- c("With interaction", "Without interaction")
AIC <- c(AIC(countmodel1), AIC(countmodel2))
resvar <- c(sum(residuals(countmodel1)^2),sum(residuals(countmodel2)^2)) #residual variance
res_sd<- c(sqrt((sum(residuals(countmodel1))^2/length(residuals(countmodel1))-2)),sqrt((sum(residuals(countmodel2))^2/length(residuals(countmodel2))-2))) # residual standard deviation
as_tibble(cbind(rownames, AIC, resvar, res_sd))
## # A tibble: 2 x 4
## rownames AIC resvar res_sd
## <chr> <chr> <chr> <chr>
## 1 With interaction 3148.44094783359 699.486557245994 0.847368170996978
## 2 Without interaction 3398.54853775285 962.462811611566 1.32245901291909
The model that includes the interaction has the lowest AIC value as well as the lowest residual variance and residual standard deviation (as compared to the model fitted with only the main effects of pas and task).
I would choose the model that includes the interaction. Partly because of the lower AIC value, residual variance and residual standard deviation, but also because it is easier to theoretically justify choosing this model. We are not really interested in predicting count from PAS rating (i.e. the participants confidence rating). Rather, we are interested in the interaction between task and pas rating, as PAS is related to the specific task at hand and we would expect PAS rating to depend on the task, considering the fact that the task variable indicates the number of potential targets that were cued (consisted of either 2,4, or 8 digit) and thus also indicates the difficulty of the task (and a more difficult task would be expected to generate a lower confidence rating).
summary(countmodel1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: count ~ pas * task + (1 + pas | subject)
## Data: data_count
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 3148.4 3232.7 -1552.2 3104.4 318
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3871 -0.7853 -0.0469 0.7550 6.5438
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 0.3324 0.5765
## pas2 0.3803 0.6167 -0.75
## pas3 1.1960 1.0936 -0.84 0.63
## pas4 2.3736 1.5407 -0.86 0.42 0.72
## Number of obs: 340, groups: subject, 29
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.03570 0.10976 36.769 < 2e-16 ***
## pas2 -0.02378 0.11963 -0.199 0.842458
## pas3 -0.51365 0.20718 -2.479 0.013166 *
## pas4 -0.77292 0.29075 -2.658 0.007853 **
## taskquadruplet 0.11490 0.03127 3.674 0.000239 ***
## tasksingles -0.23095 0.03418 -6.756 1.42e-11 ***
## pas2:taskquadruplet -0.11375 0.04605 -2.470 0.013508 *
## pas3:taskquadruplet -0.20901 0.05287 -3.954 7.70e-05 ***
## pas4:taskquadruplet -0.21500 0.05230 -4.111 3.94e-05 ***
## pas2:tasksingles 0.19536 0.04830 4.045 5.23e-05 ***
## pas3:tasksingles 0.24299 0.05369 4.526 6.02e-06 ***
## pas4:tasksingles 0.56346 0.05101 11.045 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pas2 pas3 pas4 tskqdr tsksng ps2:tskq ps3:tskq
## pas2 -0.742
## pas3 -0.829 0.613
## pas4 -0.847 0.412 0.703
## taskqudrplt -0.151 0.138 0.080 0.057
## tasksingles -0.138 0.126 0.073 0.052 0.484
## ps2:tskqdrp 0.102 -0.198 -0.054 -0.039 -0.679 -0.328
## ps3:tskqdrp 0.089 -0.082 -0.125 -0.034 -0.592 -0.286 0.402
## ps4:tskqdrp 0.090 -0.083 -0.048 -0.093 -0.598 -0.289 0.406 0.354
## ps2:tsksngl 0.098 -0.188 -0.052 -0.037 -0.342 -0.708 0.490 0.203
## ps3:tsksngl 0.088 -0.080 -0.124 -0.033 -0.308 -0.637 0.209 0.486
## ps4:tsksngl 0.092 -0.085 -0.049 -0.091 -0.324 -0.670 0.220 0.192
## ps4:tskq ps2:tsks ps3:tsks
## pas2
## pas3
## pas4
## taskqudrplt
## tasksingles
## ps2:tskqdrp
## ps3:tskqdrp
## ps4:tskqdrp
## ps2:tsksngl 0.205
## ps3:tsksngl 0.184 0.451
## ps4:tsksngl 0.507 0.474 0.427
data_count$task <- relevel(data_count$task, ref = "singles")
ggplot(data_count, aes(x = pas, y = count)) +
geom_point(aes(pas, count), color = "blue") +
facet_wrap(~ task) +
theme_bw()
bp <- ggplot(data_count, aes(x=pas, y=count, group=pas)) +
geom_boxplot(aes(fill=pas))
bp
bp + facet_grid(task ~ .)
# bp + facet_grid(pas ~ task, margins=TRUE)
A general linear mixed effect model was fitted with a poisson link function to investigate the distribution of ratings (count) as dependent on pas and task. In other words, the model predicted our outcome variable count based on the fixed effects of pas and task as well as their interaction term. For random effects, pas was modelled as a random slope for each subject, with subject modeled as random intercept. The direction of the slope (whether it is positive or negative) is highly dependent on the interaction between pas and task.
When inspecting the plots above, it seems as if the distribution of ratings across the 3 tasks (i.e, pairs, quadruplet, singles) are quite similar (with PAS 3 being the least frequent rating across all taks). Nonetheless, it seems as if the singles task had a higher count of PAS ratings 4 (i.e., “Clear experience”) which one would expect as this is the least difficult task. From just eyeballing the boxplot, it becomes more difficult to notice the differences between pairs and quadruplets tasks.
# Using the same 4 respondents as previously (3, 19, 12, 25)
data.count_four <- data_count %>%
filter(subject == "003" | subject == "019" | subject == "012" | subject == "025")
m_four <- glmer(count ~ pas * task + (pas|subject), data = data.count_four, family = poisson)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0166871 (tol = 0.002, component 1)
data.count_four %>%
ggplot() +
geom_point(aes(x = pas, y = fitted(m_four), color = "Estimated")) +
geom_point(aes(x = pas, y = count, color = "Observed")) +
facet_wrap( ~ subject)
finalmodel<- glmer(correct ~ task + (1|subject), family = "binomial", data = df_experiment)
summary(finalmodel)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ task + (1 | subject)
## Data: df_experiment
##
## AIC BIC logLik deviance df.resid
## 13636.3 13666.0 -6814.1 13628.3 12524
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7450 -1.0599 0.4791 0.6483 0.9795
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.3622 0.6018
## Number of obs: 12528, groups: subject, 29
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.11896 0.11774 9.503 < 2e-16 ***
## taskquadruplet -0.07496 0.05081 -1.475 0.14014
## tasksingles 0.16603 0.05218 3.182 0.00146 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tskqdr
## taskqudrplt -0.220
## tasksingles -0.213 0.494
Yes, task significantly explains performance (as measured by correctness) for all 3 task levels (p < .05 in all cases). In the quadruplets task, the performance is significantly worse as compared to the task pairs, but subjects’ performance in the singles task is better as compared to pairs task. This would be expected as an increase in task complexity would be expected to result in a decrease in the amount of correct answers.
finalmodel2 <- glmer(correct ~ task + pas + (1 | subject), data = df_experiment, family = 'binomial')
Task looses is significance in explaining performance. Pas seems to be a better predictor (it significantly predicts correct). This makes sense as perceptual awareness theoretically must be believed to pressupose giving the correct answer (above the threshold that would be expected by chance alone.)
finalmodel3 <- glmer(correct ~ pas + (1 | subject), data = df_experiment, family = 'binomial')
finalmodel4 <- glmer(correct ~ task*pas + (1 | subject), data = df_experiment, family = 'binomial')
rownames <- c("M1", "M2", "M3", "M4")
AIC <- c(AIC(finalmodel), AIC(finalmodel2), AIC(finalmodel3),AIC(finalmodel4))
resvar <- c(sum(residuals(finalmodel)^2),sum(residuals(finalmodel2)^2),sum(residuals(finalmodel3)^2),sum(residuals(finalmodel4)^2)) #residual variance
res_sd<- c(sqrt((sum(residuals(finalmodel))^2/length(residuals(finalmodel))-2)),sqrt((sum(residuals(finalmodel2))^2/length(residuals(finalmodel2))-2)),sqrt((sum(residuals(finalmodel3))^2/length(residuals(finalmodel3))-2)),sqrt((sum(residuals(finalmodel3))^2/length(residuals(finalmodel3))-2)),sqrt((sum(residuals(finalmodel4))^2/length(residuals(finalmodel4))-2))) # residual standard deviation
as_tibble(cbind(rownames, AIC, resvar, res_sd))
## Warning in cbind(rownames, AIC, resvar, res_sd): number of rows of result is not
## a multiple of vector length (arg 1)
## # A tibble: 5 x 4
## rownames AIC resvar res_sd
## <chr> <chr> <chr> <chr>
## 1 M1 13636.2881319358 13503.973068995 15.6567154698802
## 2 M2 12260.4130198953 12137.5692924277 13.6364434946894
## 3 M3 12256.8814750344 12138.0725449922 13.6379998755294
## 4 M4 12265.7719796063 12130.7523038036 13.6379998755294
## 5 M1 13636.2881319358 13503.973068995 13.6265684844832
anova(finalmodel, finalmodel2, finalmodel3,finalmodel4)
## Data: df_experiment
## Models:
## finalmodel: correct ~ task + (1 | subject)
## finalmodel3: correct ~ pas + (1 | subject)
## finalmodel2: correct ~ task + pas + (1 | subject)
## finalmodel4: correct ~ task * pas + (1 | subject)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## finalmodel 4 13636 13666 -6814.1 13628
## finalmodel3 5 12257 12294 -6123.4 12247 1381.4067 1 <2e-16 ***
## finalmodel2 7 12260 12312 -6123.2 12246 0.4685 2 0.7912
## finalmodel4 13 12266 12362 -6119.9 12240 6.6410 6 0.3553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1